Marching toward the eigenvalues: The Canonical Function Method and the 

Schrodinger equation 
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The Canonical Function Method (CFM) is a powerful accurate and fast method that solves the 
Schrodinger equation for the eigenvalues directly without having to evaluate the eigenfunctions. Its 
versatility allows to solve several types of problems and in this work it is applied to the solution of 
several ID potential problems, the 3D Hydrogen atom and the Morse potential. 
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I. INTRODUCTION 

The Canonical Function Method (CFM) is a powerful 
means for solving the Schrodinger equation and getting 
the eigenvalue spectrum directly in a fast and precise 
manner without computing the eigenfunctions. 

The CFM turns the two-point boundary value (TPBV) 
Schrodinger problem into an initial value problem and 
allows full and accurate determination of the spectrum. 
This is done by expressing the solution as a sum of two 
linearly independent functions (the Canonical Functions) 
with specific values at some arbitrary point belonging to 
the interval defined by the two boundaries. The integra- 
tion proceeds simultaneously from this point toward the 
left and right boundaries evaluating at each step a corre- 
sponding ratio. It stops when the difference between the 
left and right ratios is below a given desired precision. 

This work is relevant to students who have completed 
an undergraduate Quantum Mechanics course of the 
Merzbacher [l| level or graduate students whose level cor- 
responds to Landau and Lifshitz course [2| and are inter- 
ested in the eigenvalue problem of Quantum Mechanics. 

The CFM can handle a large variety of Quantum prob- 
lem problems [3J besides the eigenvalue problem making 
it an extremely versatile, fast and highly accurate. The 
evaluation of the Schrodinger operator spectrum is done 
without performing diagonalization, bypassing the evalu- 
ation of the eigenfunctions. This allows to preserve a high 
degree of numerical precision that is required in solving 
sensitive eigenvalue problems. 

It also solves the Radial Schrodinger equation over the 
infinite interval [0, oof , where singularities in the poten- 
tial at both boundaries are encountered. 

The CFM method is superior to many standard tech- 
niques that have been used to solve the Schrodinger equa- 
tion such as Numerov [f| or relaxation methods that are 
particularly tailored for solving TPBV problems (see the 
Physics Reports Review [3J). 

It is worthwhile to point out that, numerically, the 
precision gained with the bypass of intermediate diago- 
nalization operations is reminiscent of the Golub-Reinsch 
algorithm (see for instance ref. [4|) used for the singular 
value decomposition of arbitrary rectangular matrices. 



This article is organised as follows: The next sec- 
tion is a description of the CFM in ID with the appro- 
priate boundary conditions (BC). Several ID problems 
are treated: The Infinite depth potential well, the finite 
depth potential well, the Harmonic Oscillator problem, 
the Kronig-Penney potential and the double-well (sym- 
metric and asymmetric) potentials. In section III the 
CFM is applied to the 3D Schrodinger equation special- 
izing to the Radial Schrodinger equation problems and 
in particular to the Hydrogen atom and the Morse po- 
tential. Finally section IV bears our conclusions. 

In the Appendix we provide information on the differ- 
ent systems of units used in Atomic physics and quantum 
mechanics. 



II. ONE-DIMENSIONAL POTENTIAL 
PROBLEMS 

The CFM approach is based on the direct calculation 
of the eigenvalues of the Schrodinger equation defined 
over an interval [21,0:2] with a set of BC defining the 
problem: 

+ V(x)]ip(x) = Eip(x), xi<x<x 2 (1) 



2m e dx 2 



m e is the electron mass. 

Starting from a point xq £ [3:1,22] we express the 
solution as a superposition of two linearly independent 
functions a(E; 2), /3(E; 2), (the Canonical Functions) de- 
pending on the energy E and the abscissa x such that: 



y(x) = y(x )a(E; x) + y'(x )/3(E; x) 



(2) 



The CFM is based on the extraction of the eigenval- 
ues from the zeroes of the eigenvalue function F(E) de- 
fined from the saturation of the left (x — > x\) and right 
(x — >• 22) functions I -(E) and l+(E) given by the ratios 
of the canonical functions a(E;x) and /3(E;x) or their 
derivatives. 

In the general case when either y(x\) 7^ or y(x 2 ) 7^ 
we write: 



y(x) = y(x )a(E;x) +y'(x )/3(E;x) 
y'(x) = y(x )a'(E; 2) + y'(x )(3'(E; x) 



(3) 
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The canonical functions satisfy the following condi- 
tions at the starting point Xq: 

a(E; x ) = 1, a'(E; x ) = 0, 
P(E;x o )=0,p'(E;x o ) = l (4) 

Let us rewrite the system [3] at the two boundaries x = 
Xi and x = x 2 : 

V( x i) = y(xo)a(E;x 1 ) +y'(x )(3(E;x 1 ); 

y'{ Xl ) = y(x )a'(E;x 1 )+y'(x )f3'(E;x 1 ); 

y(x 2 ) = y(x )a(E;x 2 ) + y'(x o )0(E;x 2 ); 

y'(x 2 ) = y(x )a'(E;x 2 )+y'(x )(3'(E;x 2 ) (5) 

Extracting from above the left and right ratios defining 
the functions I -(E) and l+(E): 



l-(E) = 
l + (E) 



y'(xo) 



a(E; xi)y'(xx) - a'(E; xi)y(xx) 



y(xo) 
y'(xo) 



y(xo) 



p'(E;x l )y(x 1 )- 
a(E;x 2 )y'(x 2 ) 



p(E;x 1 )y'(x i y 
- a'(E;x 2 )y(x 2 ) 



P'(E;x 2 )y(x 2 ) -P(E; x 2 )y'(x 2 ) 



(6) 



In order to tackle any problem with the CFM, a num- 
ber of constraints should be explained and underlined in 
order to illustrate the methodology of getting properly 
the eigenvalue spectrum: 

• Sensitivity, stability and accuracy: 

The spectrum depends on the zeroes of F(E) = 
l+(E) — l_ (E). This subtraction might lead in some 
cases to inaccuracies because the entire spectrum 
depends on the zeroes of F(E). However, it holds 
the key of the stability of the CFM since two in- 
dependent solution sets are generated at the point 
xq , with progress inwards to the left point x\ and 
outwards toward the right point x 2 . Since both sets 
contain, in general, linear combinations of the reg- 
ular and the irregular solutions, by suitably com- 
bining them, the irregular solution is eliminated. 

• xq issue and the number of eigenvalues: 

The number of eigenvalues depend strongly on Xq. 
Thus, it should be chosen such that a tan(_E)-like 
diagram for the energy function F(E) is obtained. 
In the case we have a potential displaying a sin- 
gle minimum, xq should be close to the potential 
minimum. 

• Behaviour of the canonical functions: 

The method being sensitive to convergence of the 
marching toward the left-right boundaries X\,x 2 , 
one ought to look for similar behaviour in the 
canonical functions a(E; x) and f3(E; x) along with 
the limiting process x — > Xi and x — > x 2 since it 
controls the ratio saturation. 

• Overall aspect of the eigenvalue function: 

The eigenvalue function F(E) = l + (E) - l-(E) 
should have a regular structure of the tan(E') type, 
that is almost periodic versus \n(E). 



There are several types of BC from the eigenvalue func- 
tion defined as the difference between left and right ratio 
functions: 



F(E) = l + (E)-l_(E) = 



~y'(xo)~ 




~y'(%o)~ 


. y( x o) _ 


+ 


. y( x o) _ 



(7) 



We consider, for illustration, the following four types 
of BC: 

1. Null wavefunctions BC: 

The conditions y(xi) = y(x 2 ) = yield: 



I -(E) = lim 

X— >X 

l+(E) = lim 



a(E; x) 



x^ Xl f3(E;x)' 
a(E; x) 
J(E~xj 



(8) 



2. Null wavefunction and its derivative BC: 
The conditions y(x\) — y'(x 2 ) — yield: 



l-(E) = lim 
l + (E) 



a(E; x) 
x^txi /3(E;x)' 
a'(E;x) 



lim 

x^x2 j3'(E\x) 



3. Null derivative and the wavefunction BC: 
The conditions y'(x\) = y(x 2 ) = yield: 



(9) 



l-(E) = 
l+(E) 



a'(E-x) 
™ p'(E;xy 

= i im _4^4 

z-Hrs P(E;x) 



(10) 



4. Null derivatives BC: 

The conditions y'(xi) — y'(x 2 ) = yield: 



-(E) 



lim 



l+(E) = lim 



a'(E; x)_ 
x^x! ' P'(E;xY 
a'(E;x) 



x^x 2 P'(E;x) 



(11) 



It is remarkable that the eigenvalue function F(E) = 
1+ (E) — (E) behaves in a very peculiar way close to the 
trigonometric tan(-E) shape as displayed in Fig. [1] This 
will be explained in the next section. 



A. The Infinitely deep square well 

Let us apply the CFM to the infinite square well po- 
tential of width a defined by: 
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0.001 0.01 0.1 1 

ln(E) 

FIG. 1: Eigenvalue function versus energy displaying the first 
25 eigenvalues of the Infinite square well. The vertical lines 
indicate the eigenvalue position. The eigenvalue function has 
an approximate tan(£) shape versus the energy E. 

V(x) = 0, < x < a, V(0) = oo, V(a) = oo, meaning 
x\ = 0, X2 = a. 

The Schrodinger equation writes: 

_#iMd =mx) 0<x<a (12) 

2m e ax 

with BC: ip(0) = 0,ip(a) = 0. The exact eigenfunc- 
tions are hence given by ipn(x) = \J\ sin( 11 ^), n = 1, 2... 

yielding the exact eigenvalues as: E n = ) 2 with 

m e the electron mass. 

In order to apply, the CFM method, we first note that the 
solutions are odd or even over the interval [0, a). Working 
on the half-interval [0, |] we can start from any point xq 
and apply the general methodology albeit with a modi- 
fication regarding the matching conditions at the middle 
interval point. 

In the odd-mode case (null wavefunctions at both 
boundaries x\ = 0, x-i = |): 



Index 


CFM 


Exact 


1 


1.6006952(-3) 1.6000001(-3) 


3 


1.4406255(-2) 1.4400001(-2) 


5 


4.0017359(-2) 4.0000003(-2) 


7 


7.8434058(-2) 7.8400001(-2) 


9 


0.1296562 


0.1296000 


11 


0.1936839 


0.1936000 


13 


0.2705172 


0.2704000 


15 


0.3601561 


0.3600000 


17 


0.4626004 


0.4624000 


19 


0.5778504 


0.5776000 


21 


0.7059059 


0.7056000 


23 


0.8467670 


0.8464000 


25 


1.000433 


1.000000 


2 


6.4027691(-3) 6.4000003(-3) 


4 


2.5611134(-2) 2.5600001(-2) 


6 


5.7624962(-2) 5.7600003(-2) 


8 


0.1024444 


0.1024000 


10 


0.1600693 


0.1600000 


12 


0.2304998 


0.2304000 


14 


0.3137361 


0.3136000 


16 


0.4097775 


0.4096000 


18 


0.5186247 


0.5184000 


20 


0.6402774 


0.6400000 


22 


0.7747357 


0.7744000 


24 


0.9219995 


0.9216000 



TABLE I: First twenty-five odd and even quantum levels of 
the infinite square well potential given by the CFM along with 
exact results. The well width a is chosen in a way such that 
the eigenvalue is 1 when the level index is 25. The numbers 
in parenthesis represent a power of 10. All eigenvalues are in 
Atomic units (see Appendix). 



B. The Finite depth square well 



at(E; x) 
""o ' p(E; x) ' 
a(E; x) 



I- (E) = lim 

l+{E)=Um 

*->•§ p(E;x) 

F {E)=l+{E)-l_{E) 



(13) 



whereas in the even-mode case we have (null wave- 
function at left boundary x\ = 0, null derivative at right 
boundary X2 — §): 



a(E; x) 
x^h P{E;xY 
a'{E;x) 



I -(E) = lim - 
l+(E)= lim 

an^f P'(E;x) 

Fe(E) = l+(E) — l-(E) 



(14) 

Solving successively F a (E) and F e (E) for the odd 
modes and the even modes, we get Table HI 



The finite depth square well potential of width a is 
defined by: V(x) = 0, for < x < a; V(x) = 
Vq] for x > a or x < 0. 

As in the Infinite depth case, the potential is sym- 
metric with respect to the well center |, implying 
that we have odd and even modes. Therefore we take 
x\ = —oo,X2 = % which means that we march to the left 
through the potential step until we observe the nulling 
of the wavefunction, whereas the marching to the right 
results at half the potential well width a in odd or even 
modes. More explicitly: 

In the odd-mode case (null wavefunctions at both 
boundaries x\ = — oo,:^ = §): 



1 _ {E )= lim -^4; 
v ; z^-°o P(E-xY 
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1+{E) = lim - 



a(E; x) 



F (E) = 1+{E) - l_(E) 



(15) 



whereas in the even-mode case we have (null wavefunc- 
tion at left boundary x% = — oo, null derivative at right 
boundary X2 — §): 



1-{E)= lim 
= lim 



a(-E; x) 
v "/3(i?;x) ; 



0'(E;xY 
F e (E) = l+(E)-l_(E) 



(16) 



Solving successively F a (E) and F e (E) for the odd 
modes and the even modes, we get the following table [TT1 



Index 


CFM 


Exact 


1 


1.6482281(-3) 1.6475233(-3) 


3 


1.4832322(-2) 1.4826014(-2) 


5 


4.1191306(-2) 4.1173782(-2) 


7 


8.0705732(-2) 8.0671579(-2) 


Q 


0.1333447 


0.1332882 


11 


0.1990615 


0.1989772 


13 


0.2777881 


0.2776708 


15 


0.3694246 


0.3692691 


17 


0.4738183 


0.4736193 


19 


0.5907167 


0.5904695 


21 


0.7196453 


0.7193471 


23 


0.8594448 


0.8590948 


2 


6.5926472(-3) 6.5898113(-3) 


4 


2.6365897(-2) 2.6354689(-2) 


6 


5.9305709(-2) 5.9280563(-2) 


8 


0.1053872 


0.1053425 


10 


0.1645720 


0.1645023 


12 


0.2368039 


0.2367037 


14 


0.3220005 


0.3218648 


16 


0.4200394 


0.4198628 


18 


0.5307267 


0.5305039 


20 


0.6537226 


0.6534503 


22 


0.7883219 


0.7879974 


24 


0.9322464 


0.9318770 



TABLE II: First 24 odd and even quantum levels of the finite 
depth square well potential given by the CFM along with 
exact results. The numbers in parenthesis represent a power 
of 10. The barrier height Vo = 1 and all eigenvalues are in 
Atomic units (see Appendix). 



fife 

The exact eigenvalues E n — , drawn from Landau- 
Lifshitz book [2j are given by the solutions k n (n — 
1,2...) of the transcendental equation: 



sin 



flkr 



= —(mr — k n a): < k n < 



2mVb 



(17) 




ln(E) 



FIG. 2: Eigenvalue function versus energy displaying the first 
15 eigenvalues of the finite square well. The vertical lines 
indicate the eigenvalue position. The eigenvalue function has 
an approximate tan(_E) shape versus the energy E. 



The number of levels gives us a the well width in the 
following way: since the sin" 1 term is bounded by ^, the 

= 1 + 



largest level index n max is given by n ri 



2m Vp 



hence a 



7rh(n r) 



-1) 



v / 2mVb 



C. The harmonic oscillator 

The harmonic oscillator potential defined by: V(x) — 
■^kx 2 is symmetric with respect to the origin x — im- 
plying as before that the solutions are given by odd and 
even parity modes. The boundaries for this problem are: 

Xi = — OO, X2 = 0. 

In the odd-mode case (null wavefunctions at both 
boundaries x\ — — oo, X2 = 0): 



I -(E) = lim - 



l + (E) - 
Fo(E) 



a(E; x) 
J(E~x) ] 
a(E; x) 
x^o P(E;x)' 
l + (E)-l_(E) 



lim 



(18) 



whereas in the even-mode case we have (null wavefunc- 
tion at left boundary x\ = — oo, null derivative at right 
boundary X2 = 0): 



l_(E) = lim - 

x — y — oc 



a(E; x) 
J(E^xj ] 

l + (E) = lim -f E > X l 
+ V ' x^o P'(E; X y 

F e (E) = l+(E)-l_(E) 



(19) 



Solving successively F Q (E) and F e (E) for the odd 
modes and the even modes, we get the following table Mil 

The exact eigenvalues E n — hoj n (n + |) Q allow us to 
pick the energy of the highest level as 1 (in Atomic units) 
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Index 


CFM 


Exact 


1 


5.8842082(-2) 5.8823533(-2) 


3 


0.1372950 


0.1372549 


5 


0.2157480 


0.2156863 


7 


0.2941877 


0.2941177 


9 


0.3726121 


0.3725490 


11 


0.4510336 


0.4509804 


13 


0.5294515 


0.5294118 


15 


0.6078746 


0.6078432 


17 


0.6863770 


0.6862745 


19 


0.7654451 


0.7647059 


21 


0.8468009 


0.8431373 


23 


0.9335056 


0.9215686 


25 


1.028182 


1.000000 





1.9617772(-2) 1.9607844(-2) 


2 


9.8067097(-2) 9.8039217(-2) 


4 


0.1765225 


0.1764706 


6 


0.2549707 


0.2549020 


8 


0.3334008 


0.3333333 


10 


0.4118234 


0.4117647 


12 


0.4902429 


0.4901961 


14 


0.5686612 


0.5686275 


16 


0.6471026 


0.6470588 


18 


0.7257724 


0.7254902 


20 


0.8056628 


0.8039216 


22 


0.8892854 


0.8823529 


24 


0.9797482 


0.9607844 



TABLE III: Ground state (zero index) and first twenty-five 
odd and even excited states of the harmonic oscillator poten- 
tial given by the CFM along with exact results. The numbers 
in parenthesis represent a power of 10. The oscillator elastic 
constant was chosen such that level 25 had value 1 in Atomic 
units. All eigenvalues are in Atomic units (see Appendix). 



from which we select the value of ujq = 




and hence 



the elastic constant k. 



D. Periodic potential: The Kronig-Penney problem 

The Kronig-Penney potential is often used in the de- 
scription of the electronic properties of crystals. It is 
based on a piecewise constant potential (see fig. QJ for 
which we can apply the same methodology of marching 
to the left and to the right and comparing correspond- 
ing ratios in order to get the eigenvalues. The latter are 
now dispersive which means they depend on a wavevec- 
tor reflecting the translational symmetry of the system 
(Bloch theorem). The CFM must be extended to the 
complex case since previously all the wavefunctions we 
use and derive were real. It is straightforward to extend 
the marching method as well to the complex wavefunc- 
tion case as explained below. 

v, v. 

E 



-b a x 

FIG. 4: Periodic piecewise constant potential V(x) displaying 
alternating regions of V — and V = Vo with periodicity a+b. 
The energy bands are obtained for E < Vo- In the case we 
let Vo — > oo and b — > the barriers become delta functions 
sitting on a periodic lattice with parameter a. 

Defining a unitcell with extreme boundaries —b and 
a we write the general CFM definitions in the complex 
case: 




FIG. 3: Eigenvalue function versus energy displaying the first 
26 eigenvalues of the harmonic oscillator potential. The ver- 
tical lines indicate the eigenvalue position. The eigenvalue 
function has an approximate tan(E') shape versus the energy 
E. 



y(-b) = y(xo)a(E;-b)+y'(x )l3(E;-b) 

y'(-b) = y(x )a'(E;-b) + y'(x )f3'(E;-b) 

y(a) = y(x )a(E;a) + y'(x )/3(E;a) 

y\a) = y(x )a'{E;a)+y'{x )P'(E;a) (20) 

The energy E is considered as smaller than Vq. Using 
Bloch theorem [||, in the above equations: 

y(a) = y{—b) exp[ik(a + b)], y'(a) = y'{—b) exp[ik(a + b)], 

(21) 

we get the complex ratio functions: 



1-{E) = 
l + (E) = 



y'(xo) 



_ ja(E; -b) - a{E;a) _ 
" [}{E-a)- 1 l3{E--by 
-ya'(E;-b) - a'(E;a) 
P'{E-a)- 1 (3'{E;-by 
F{E) = l + (E) - /_(£) 



(22) 



where 7 = exp[ifc(a + b)]. This yields the dispersion 
relation for the energy eigenvalue e„(fc) with n the band 
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index: 

[ya(E; -b) - a(E; a)][p'(E; a) - j/3'(E; -b)] - 
[/3(E; a) - 7/3(£; -b)] [ja'{E; -6) - a'(E; a)} = 0(23) 

We compare the above to the standard dispersion re- 
lation [f|: 

/02 ,„2 

• smh(Qb) sin(fta) +cosh(<3&) cos na — cos k(a + b) 



2Qk 



(24) 

where Q is defined as the (pure imaginary) wavevec- 
tor inside the potential barrier. Recall that the en- 
ergy E < Vq hence the wavefunction within the bar- 
rier is of the form exp(±Qx), i.e. when V(x) = Vq, 

ji^n 2 r 2 2 

Vq = o + whereas k is the real free wavevector 

u 2m e 2m e 

outside the barrier i.e. when V(x) = the wavefunction 
is of the form exp(±inx)). 

When we let Vq — > oo and b — » such that Vo& remains 
finite, the piecewise constant potential V(x) is trans- 
formed into a periodic array of S functions gS(x — na) 
with lattice parameter a. g is the strength of the delta 
function potential and n € Z a relative integer. 

We can formally write the potential as: 



n — +oo 



(25) 



and consider a single interval extending over the unit cell 
with boundaries x% = and £2 = a. Since at the left 
boundary X\ = we have a i5 function potential, stan- 
dard quantum mechanics [l| tell us that the wavefunction 
derivative y'{x) jumps across the S function potential, 
such that: 



Bloch theorem [6j] transforms this equation into: 

y'(0+) - y'(a-)exp(-ika) = g—^y(0) 

n 

The left ratio (complex) is thus obtained as: 

y'jxo) _ 
y{xo) 



(26) 



(27) 



l-(E) = 



-a'(E; 0) + a'(E; a) exp(-ika) + g^a{E; 0) 
P'{E; 0) - P'{E; a) exp(-ifca) - g^fi(E; 0) 



(28) 



The right ratio is determined from Bloch theorem link- 
ing the right boundary xi = a to the left boundary 
x\ = 0: y(a) = y(0) exp(ifca): 

7 = ^( x o) = Q (^i a ) -exp(zfca)a(£;0) 
+ l ' y{x ) exp(ika)p{E;0) - 0(E; a) { > 



The dispersion relation is obtained as before from the 
zeroes of: 



F(£) = *+(£)- /_(£) = 



a(E; a) — exp(ika)a(E; 0) 



9^ 



exp(ika)p{E;0) -0(E;a) 
a'{E; 0) + a'(E; a) exp(-ika)a(E; 0) 



(30) 



P'(E; 0) - P'(E; a) exp(-ika) - g^f3(E; 0) 

Indeed, the dispersion relation [6] obtained from the 
limiting process Vq — > 00 and b — > is [y] : 



Q 2 b . 

— — sm(Ka) + cos na 
2k 



cos ka 



(31) 



can be straightforwardly obtained from the deriva- 
tive jump condition (eq. |27|) and Bloch theorem y(a) = 
y(0) exp(ika). Starting with the wave function y(x) — 
Aexp(inx) + B exp(— inx), defined over the unit cell 
x G]0, a[ and using both aforementioned conditions yields 
the dispersion relation: 



■ sin(«;a) + cos kg = cos ka 



(32) 



Comparing both dispersion relations yields finally the 

value of the strength of the 5 function potential as g = 

Q 2 bh 2 



(a) 



(b) 



(c) 



FIG. 5: (Color on line) Exact bands (in green or " x") for the 
Kronig-Penney model and comparison with the CFM bands 
(in red or " +" ) obtained from the dispersion relation obtained 
from ea. l30l and l28l (a) is for a single band, (b) and (c) are for 
3 and 5 bands respectively with a — 2.22, 6.66, 11.12 Atomic 
units. In all cases, the value of xq is 1 and the strength of the 
potential g = 1 in Atomic units. 



In figure [5j exact bands are compared to the CFM 
bands. The lattice parameter in each case is deter- 
mined from the number of bands ns we want to cal- 
culate according to the formula: a = since the 

largest wavenumber is k rnax — and we select the 



largest energy E max — as 1. This is why we use 

a = 2.22,6.66,11.12 Atomic units for the ub — 1,3,5 
respective number of bands. It is remarkable to observe 
how the CFM results lie exactly on top of the exact re- 
sults. The starting value Xq is chosen in a way such 
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that we get the right number of bands riband- Spurious 
bands might appear due to a bad starting value Xq be- 
cause the nature of the CFM dispersion relationl30l differs 
with respect to the dispersion relations eq. [3T] and eg I3"2"1 
The latter eqj32] allow the exact determination of the free 
wavevector k from a given Bloch wavevector k and the 
exact band energy is obtained from E — tj- 2 -. In sharp 
contrast, the CFM dispersion relation [30l yields directly 
the band energy without going through the determina- 
tion of an intermediate wavevector k. 



E. Double well potential over an infinite interval 

Double-minimum Potential Well (DPW) problems de- 
fined over the semi-infinite interval [0, oof are interesting 
to solve as they arise in many areas of Atomic, Molecular 
and even Solid State physics. When two-dimensional 
electron layers (such as in heterostructures involving 
semiconductors) are placed in perpendicular electric and 
magnetic fields, a potential well with two minima, for 
electronic motion normal to the surface, arises. 

A DPW can be symmetric or asymmetric and one has 
to adapt in each case the appropriate BC imposed by the 
CFM. 

A symmetric DPW is the Double Gaussian potential 
investigated by Hamilton and Light [7] given by: 

V(x) = -D[exp(-Q(x - r a f) + exp(-ft(z + r a ?)] 

The values of the parameters: D, ft, r a are respectively: 
12.0,0.1,5.0 in standard atomic units (see Appendix) such 
that h = 1, m e = 1. 

An elaborate method used by Hamilton and Light Q 
based on Distributed Gaussian Basis sets borrowed from 
Quantum Chemistry Techniques gives the eigenvalues 
listed in tabic [TV] The CFM results for all the 24 levels 
in table IIVI proves once again that it is able to find all 
levels with speed and accuracy from a simple marching 
approach. 

The asymmetric DWP introduced by Johnson consists 
of the sum of a Morse (see next section) and a Gaussian 
potentials such that: 

V(x) = D[l - exp(-B(x - r a ))] 2 + A exp (— C(x — r ) ) 

The values of the parameters A, B, C, D 7 r a ,rb are (fol- 
lowing Johnson Q) in (cm -1 , A system of units) are: 10 4 
cm" 1 , 1.54 A" 1 , 200.0 A" 2 , 31250.0 crn^ 1 , 1.5 A, 1.6 A 
respectively. 

Eigenvalues for the asymmetric double minimum po- 
tential problem are given in table |V] and a comparison 
between Johnson's ;8!] results and the CFM are displayed 
below. 



III. THE CANONICAL FUNCTION METHOD 
AND THE 3D RADIAL SCHRODINGER 
EQUATION 



After having discussed the CFM in the ID case, we 
move on to the treatment of the Radial Schrddinger 
Equation (RSE) . The mathematical difficulty of the RSE 
lies in the fact it is a singular boundary value problem 
(SBVP). problem. This stems from the boundary 
conditions over the infinite interval [0,oo[ , with the 
double requirement of regularity near the origin (r ~ 0) 
where the potential is large and near infinity (r — > oo) 
where the potential is very small. The CFM turns 
it into a regular initial value problem and allows the 
full determination of the spectrum of the Schrddinger 
operator bypassing the evaluation of the eigenfunctions. 

The partial wave form of the RSE is written as: 



Index 


CFM 




Hamilton and Light 


1 


-11.250 421 409 


-11.245 199 313 


3 


-9.779 


225 


834 


-9.773 496 


902 


5 


-8.387 


719 


137 


-8.381 307 


491 


7 


-7.079 


412 


929 


-7.072 038 


846 


9 


-5.858 


811 


221 


-5.849 940 





11 


-4.732 


171 


001 


-4.720 509 


6 


13 


-3.709 


113 


861 


-3.690 475 


6 


15 


-2.801 


628 


760 


-2.763 219 


7 


17 


-2.000 


566 


637 


-1.924 577 




19 


-1.255 


332 


005 


-1.149 254 




21 


-0.561 


216 


170 


-0.457 88 




23 


-0.045 


810 


537 


-0.003 41 







-11.250 418 469 


-11.245 199 313 


2 


-9.779 


202 


594 


-9.773 496 


902 


4 


-8.387 


701 


732 


-8.381 307 


510 


6 


-7.079 


415 


041 


-7.072 039 


562 


8 


-5.858 


805 


474 


-5.849 958 


02 


10 


-4.732 


231 


858 


-4.720 829 


36 


12 


-3.709 


907 


559 


-3.694 518 


38 


14 


-2.807 


436 


691 


-2.798 251 


92 


16 


-2.022 


064 


904 


-2.089 661 


3 


18 


-1.293 


090 


067 


-1.462 202 


9 


20 


-0.601 


483 


056 


-0.771 081 




22 


-0.067 


153 


689 


-0.177 181 





TABLE IV: Computed odd and even eigenvalues for the sym- 
metric double Gaussian well potential. The numbers at left 
are the levels computed with the CFM; on the right the levels 
obtained by Hamilton and Light 0]. Note the deterioration 
of accuracy of the Hamilton and Light results as the index 
increases because of the approach of the continuum. 
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Index 


Johnson 


CFM 





1302.500 


1302.498 972 


1 


3205.307 


3205.303 782 


2 


4227.339 


4227.336 543 


3 


5144.251 


5144.243 754 


4 


6064.241 


6064.225 881 


5 


7092.679 


7092.664 815 


6 


7614.622 


7614.603 506 


7 


8911.545 


8911.513 342 


8 


9095.696 


9095.679 497 


9 


10208.350 


10208.318 142 


10 


10869.289 


10869.255 077 


11 


11482.475 


11482.457 956 


12 


12353.799 


12353.766 422 


13 


12972.473 


12972.453 117 


14 


13690.455 


13690.436 602 


15 


14435.350 


14435.321 044 



TABLE V: Eigenvalues in cm 1 of the Johnson asymmetric 
DWP consisting of the sum of a Morse and a Gaussian poten- 
tials V(r) = D[l - exp(-B(x - r a ))] 2 + Aexp (-C(x - r b f). 
with A= 10 4 cm -1 , B= 1.54 A -1 , G= 200.0 A -2 , D= 31250.0 
cm -1 , and r a = 1.5 A, rj,= 1.6 A. Johnson Q] results are com- 
pared to the CFM. 



h 2 d 2 Ri(E;r) 
2fi dr 2 



V(r) 



h 2 + 



2fi 



Ri{E;r) = 
ERi(E;r\33) 



where /j, is the reduced mass and Ri (E; r) is the reduced 
probability amplitude for orbital angular momentum I 
and eigenvalue E. 

The BC are: 



lim Ri(r) = 0; lim Ri(r) 

r—>0 t— >-+oo 







(34) 



The CFM consists of writing the general solution y(r) 
representing the probability amplitude Ri(E; r) as a func- 
tion of the radial distance r in terms of two linearly in- 
dependent basis functions a(E; r) and P(E; r) for some 
energy E. 

Generally, the RSE is rewritten in a system of units such 
that h = 1, 2/x = 1 (see Appendix on units): 



d 2 y(r) 
dr 2 



V(r) 



1(1 + 1) 



- E 



y(r) (35) 



At a selected distance ro, a well defined set of initial 
conditions are satisfied by the canonical functions and 
their derivatives ie: a(E;ro) = 1 with a'(E;ro) = and 
P(E\ ro) = with /3'(E; r ) = 1. Thus we write as in the 
ID case: 



y(r) = y(r )a(E; r) + y'(r )(3(E; r) 



(36) 



The method of solving the RSE is to proceed from ro 
simultaneously towards the origin (r — > 0) and towards 
infinity (r — > oo). 

When the integration is performed, the ratio of the r 
dependent canonical functions is monitored until satura- 
tion with respect to r is reached at both limits (r — > 
and r — » oo). The saturation of the g7gpj ratio with r 
yields a position independent eigenvalue function F(E): 



F(E) = l + (E)-l_(E) 



|Y(ro)~ 




|Y(ro)l 


. y( r o) _ 


+ 


. y( r o) . 



(37) 



The tan(-E) shape of F(E) provides a deep insight into 
the physical significance of the CFM method. The lat- 
ter transforms a SB VP from the open interval [0,oo[ to 
a finite interval [rj e /t f ^Hgftt] defined by the saturation 
coordinates of the ratio functions. This means the CFM 
maps an arbitrary potential V(r) onto the infinite square 
well problem in the finite interval \riefti fright] resulting 
in an eigenvalue function F(E) with a tan(-E) pattern as 
we saw in Section II (see also ref. [1]). 



A. The Hydrogen atom spectrum 

The Coulomb potential is a crucial case to test the ac- 
curacy and reliability of the CFM given by the Hydrogen 
atom problem. 

The CFM results are shown in Table. I VII along with 
the exact analytical results and it is remarkable to notice 
that all digits (calculated by CFM and analytically) are 
rigorously same. 



B. The Morse potential 

The classical Morse potential is the simplest model 
for the evaluation of cell vibrational spectra of diatomic 
molecules. 

The Morse potential is given by: 



V(r) = D[l - cxp(-a{r - r e })Y 



D 



(38) 



with the values D, a, r e equal respectively to 188.4355, 
0.711248, 1.9975 in atomic units (see Appendix). The 
analytic expression for the levels is: 



E n 



a 2 h 2 q/2^D 




(39) 



Hence the number of levels 



with max n < 

is given by: N = 

Working with units such that h = 1 and 2/i = 1, the 
Morse potential and the eigenvalue function F(E) are 
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Index CFM (Ry) Exact (Ry) 



1 




-1.00000 




-1.00000 


2 


- 


0.250000 


- 


0.250000 


3 


- 


0.111111 


- 


0.111111 


4 


-6 


25000 


-2) 


-6 


25000 


-2) 


5 


-4 


00000 


-2) 


-4 


00000 


-2) 


6 


-2 


77778 


-2) 


-2 


77778 


-2) 


7 


-2 


04082 


-2) 


-2 


04082 


-2) 


8 


-1 


56250 


-2) 


-1 


56250 


-2) 


9 


-1 


23457 


-2) 


-1 


23457 


-2) 


10 


-1 


00000 


-2) 


-1 


00000 


-2) 


11 


-8 


26446 


-3) 


-8 


26446 


-3) 


12 


-6 


94444 


-3) 


-6 


94444 


-3) 


13 


-5 


91716 


-3) 


-5 


91716 


-3) 


14 


-5 


10204 


-3) 


-5 


10204 


-3) 


15 


-4 


44445 


-3) 


-4 


44445 


-3) 


16 


-3 


90625 


-3) 


-3 


90625 


-3) 


17 


-3 


46021 


-3) 


-3 


46021 


-3) 


18 


-3 


08642 


-3) 


-3 


08642 


-3) 


19 


-2 


77008 


-3) 


-2 


77008 


-3) 


20 


-2 


50000 


-3) 


-2 


50000 


-3) 


21 


-2 


26757 


-3) 


-2 


26757 


-3) 


22 


-2 


06612 


-3) 


-2 


06612 


-3) 


23 


-1 


89036 


-3) 


-1 


89036 


-3) 


24 


-1 


73611 


-3) 


-1 


73611 


-3) 



TABLE VI: Energy levels of the Hydrogen atom. Middle col- 
umn values are the CFM results whereas the last column val- 
ues are the corresponding exact analytically obtained values. 
The numbers in parenthesis represent a power of 10. 

displayed in figJS] and fig. [7] respectively. Table IVIII con- 
tains the levels calculated by CFM and compared to the 
analytical analytical results. As in all previous cases, the 
agreement is perfect and the full set of levels (N = 19) 
are found as predicted analytically. 



200 r- 
150 - 
100 - 




-20(1 



2 4 6 8 10 



FIG. 6: Morse potential V{r) = D[l - exp(-a{r - r e })] 2 ~D 
with parameters D = 188.4355, a = 0.711248, r e = 1.9975. 



Index 


CFM 


Exact 


1 


-178.798248 


-178.798538 


2 


-160.282181 


-160.283432 


3 


-142.778412 


-142.78006 


4 


-126.287987 


-126.288445 


5 


-110.807388 


-110.808578 


6 


-96.3395233 


-96.3404541 


7 


-82.8832169 


-82.884079 


8 


-70.4389801 


-70.4394531 


9 


-59.0056 


-59.0065727 


10 


-48.5851288 


-48.5854378 


11 


-39.1754532 


-39.1760521 


12 


-30.77771 


-30.7784157 


13 


-23.3919983 


-23.3925247 


14 


-17.0183048 


-17.018383 


15 


-11.6557436 


-11.6559868 


16 


-7.3050122 


-7.30533791 


17 


-3.9661877 


-3.9664371 


18 


-1.6390723 


-1.63928342 


19 


-0.3238727 


-0.32387724 



TABLE VII: Energy levels of the Morse potential V(r) = 
D[l — exp(— a{r — r e })] 2 — D with parameters D — 
188.4355, a = 0.711248, r e = 1.9975. Middle column values 
are the CFM results whereas the last column values are the 
corresponding exact analytically obtained values. Units are 
such that h = 1 and 2/i = 1. 



0.1 1 10 100 

ln(IEI) 

FIG. 7: Behavior of the eigenvalue function F(E) with abso- 
lute value of energy on a semi-log scale for the Morse poten- 
tial. 



IV. CONCLUSIONS 

The CFM is a very powerful, fast and accurate method 
that is able to evaluate the eigenvalue spectrum without 
having to determine first or simultaneously the eigenfunc- 
tions. 
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The CFM has been tested succesfully in a variety of 
potentials [|| and gives accurate results for bound and 
free states. The tunable accuracy of our method allows 
to evaluate eigenvalues close to the ground state as well 
as close to highly excited states near the continuum limit 
to a large number of digits without any extrapolation. 



, which has dimensions ML 3 T~ 2 , and another quantity 
that appears all over in quantum physics is Ti which has 
dimensions ML 2 T _1 ; so it is convenient to choose units 
of length and time such that 47reo = 1 and Ti = 1. 

Using dimensional analysis, the atomic unit of length 

is: 



The CFM compares favorably with many different 
elaborate techniques based on expansion over basis func- 
tions (such as Gaussian Q, Quantum Chemistry in- 
spired basis functions [3]) or functional expansions (Nu- 
merov @, High-order Taylor pj...). The CFM approach 
remains the same despite the wide variability of the men- 
tioned problems. 

The CFM method used gives the right number of all 
the levels and the variation of the eigenvalue function 
F(E) definitely determines the total number of levels. 
Generally in order to avoid potential singularities, Taylor 
series expansion are made to a given order dictated by 
the required accuracy (as described in ref [9J). 

Since the CFM bypasses the calculation of the eigen- 
functions, it avoids losing accuracy associated with the 
numerical calculation specially with rapidly oscillating 
wave functions of highly excited states. This is specially 
needed in the study of the sensitive problem of Rydberg 
states in Atomic physics or the determination of vibra- 
tional spectra of cold (weakly-bound) molecules... 
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APPENDIX 
Atomic and other units 

In atomic and molecular physics, it is convenient to 
use the elementary charge e, as the unit of charge, and 
the electron mass m e as the unit of mass (despite the 
fact that in some cases the proton mass, m p , or the 
unified mass unit amu, is more convenient). Electrostatic 
forces and energies in atoms are proportional to e 2 /47reo 



a B 



m e (e 2 / 4ire ) '' 



(40) 



called the Bohr radius, or simply the bohr (0.529 A), 
because in the "Bohr model" the radius of the smallest 





J 


eV 


Hz 


; m, p i 

-1 

cm 


J 


1 


6.24151. 10 18 


1.50919.10 33 


5.03411. 10 22 


eV 


1.60219. HT 19 


1 


2.41797.10 14 


8.06547.10 3 


Hz 


6.62619. HP 34 


4.13570.1CT 15 


1 


3.33564.HT 11 


cm -1 


1.96648. HT 23 


1.23935.10 ~ 4 


2.99792.10 10 


1 



TABLE VIII: Conversion table for the energy expressed in J, 
eV, Hz and cm -1 . In MKS the Joule is preferred whereas 
physicists in general use eV or Hz. Spectroscopists and 
Chemists use rather the cm -1 . 



In full quantum theory the particles do not follow an or- 
bit but possess wavefunctions and the expectation value 
of the electron-proton distance in the Hydrogen ground 
state is exactly (1 + ^)«s- 

The atomic unit of energy is the Hartree (27.2 eV) 
given by: 



E h 



1 



(41) 



The unit of time is ti/ E^. 

The Hartree is twice the ground state energy of the 
Hydrogen atom hO- + Eh equal to the Rydberg 

(13.6 eV). In atomic and molecular spectroscopy, one uses 
rather the cm -1 an energy corresponding to a wavelength 
of 1cm or sometimes a frequency unit, the Hz. We refer 
the reader to table IVIHl where conversion factors between 
the different energies are given. 
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